import numpy as np
import array
import pandas as pd
import matplotlib.pyplot as plt
import data_entry2

# DMM measured resistance
R2 = 9970 # Ohms #99.99 kOhms, so resolution is 0.01kOhms, accuracy is 0.5%+3*unit
dR2 = 9970*0.005+3*(0.01*1000)

R1_start = 2530 # Ohms #9.999 kOhms, so resolution is 0.001 kOhms, accuracy is 0.5%+3*unit
dR1_start = 2530*0.005+3*(0.001*1000)

R3 = 1.004e6 # Ohms #9.999 MOhms, so resolution is 0.001 MOhms, accuracy is 1.5%+3*unit
dR3 = 1.004e6*0.015 + 3*(0.001*1e6)

# You may notice that i use the value of some uncertainty directly, because i want to hold all values in one counting method
print(f"DMM: 10 fixed kOhms(R2) resistor is {R2} +/- 80 Ohms, 10 kOhms variable resistor(R1_start) is {R1_start} +/- 20 Ohms, 1 MOhms resistor(R3) is 1000000 +/- 20000 Ohms")

# DMM measured capacitance
C1 = 9.576 # nF #9.999 nF, so resolution is 0.001 nF, accuracy is 5%+20*unit
dC1 = 9.576*0.05+20*(0.001)

C2 = 206.2 # F #999.9 uF, so resolution is 0.1 uF, accuracy is 2%+5*unit
dC2 = 206.2*0.02+5*(0.1)
print(f"DMM: 0.01uF Capacitor(C1) is {C1:.1f} +/- {dC1:.1g} nF, 220uF Capacitor(C2) is {C2:.0f} +/- {dC2:.1g} uF")

DMM: 10 fixed kOhms(R2) resistor is 9970 +/- 80 Ohms, 10 kOhms variable resistor(R1_start) is 2530 +/- 20 Ohms, 1 MOhms resistor(R3) is 1000000 +/- 20000 Ohms
DMM: 0.01uF Capacitor(C1) is 9.6 +/- 0.5 nF, 220uF Capacitor(C2) is 206 +/- 5 uF

# Amp = 200 mV with Voffset = 0 
Vin1 = 0.204 # V
uVin1 = (0.006)/2 # V
Vout1 = 0.960 # V
uVout1 = (0.980 - 0.940)/2 # V

# Amp = 300 mV with Voffset = 0 
Vin2 = 0.306 # V
uVin2 = (0.006)/2 # V
Vout2 = 1.46 # V
uVout2 = (1.48 - 1.44)/2 # V

# Amp = 400 mV with Voffset = 0 
Vin3 = 0.400 # V
uVin3 = (0.006)/2 # V
Vout3 = 1.97 # V
uVout3 = (1.99 - 1.96)/2 # V

# Amp = 500 mV with Voffset = 0 
Vin4 = 0.500 # V
uVin4 = (0.006)/2 # V
Vout4 = 2.50 # V
uVout4 = (2.50-2.47)/2 # V

# Amp = 600 mV with Voffset = 0 
Vin5 = 0.598 # V
uVin5 = (0.006)/2 # V
Vout5 = 2.96 # V
uVout5 = (2.97 - 2.95)/2 # V

# Amp = 700 mV with Voffset = 0 
Vin6 = 0.700 # V
uVin6 = (0.006)/2 # V
Vout6 = 3.47 # V
uVout6 = (3.48 - 3.46)/2 # V

# Amp = 800 mV with Voffset = 0 
Vin7 = 0.796 # V
uVin7 = (0.006)/2 # V
Vout7 = 3.97 # V
uVout7 = (3.98-3.96)/2 # V

# Amp = 900 mV with Voffset = 0 
Vin8 = 0.892 # V
uVin8 = (0.006)/2 # V
Vout8 = 4.43 # V
uVout8 = (4.44-4.42)/2 # V

# Amp = 1000 mV with Voffset = 0 
Vin9 = 0.992 # V
uVin9 = (0.006)/2 # V
Vout9 = 4.92 # V
uVout9 = (4.94-4.90)/2 # V

# Amp = 1100 mV with Voffset = 0 
Vin10 = 1.09 # V
uVin10 = (0.006)/2 # V
Vout10 = 5.44 # V
uVout10 = (5.46-5.42)/2 # V

# Amp = 1200 mV with Voffset = 0 
Vin11 = 1.19 # V
uVin11 = (0.006)/2 # V
Vout11 = 5.92 # V
uVout11 = (5.94-5.90)/2 # V

# Amp = 1340 mV with Voffset = 0 
Vin12 = 1.33 # V
uVin12 = (0.006)/2 # V
Vout12 = 6.56 # V
uVout12 = (0.080)/2 # V

# Amp = 1440 mV with Voffset = 0 
Vin13 = 1.43 # V
uVin13 = (0.006)/2 # V
Vout13 = 6.92 # V
uVout13 = (0.080)/2 # V

# Amp = 1540 mV with Voffset = 0
Vin14 = 1.52 # V
uVin14 = (0.006)/2
Vout14 = 7.24 # V
uVout14 = (0.080)/2

V1 = [Vin1, Vin2, Vin3, Vin4, Vin5, Vin6, Vin7, Vin8, Vin9, Vin10, Vin11, Vin12, Vin13, Vin14]
V2 = [Vout1, Vout2, Vout3, Vout4, Vout5, Vout6, Vout7, Vout8, Vout9, Vout10, Vout11, Vout12, Vout13, Vout14]
uV1 = [uVin1, uVin2, uVin3, uVin4, uVin5, uVin6, uVin7, uVin8, uVin9, uVin10, uVin11, uVin12, uVin13, uVin14]
uV2 = [uVout1, uVout2, uVout3, uVout4, uVout5, uVout6, uVout7, uVout8, uVout9, uVout10, uVout11, uVout12, uVout13, uVout14]
plt.errorbar(V1,V2, yerr=uV2, marker=".", linestyle="",label="Measured data")
plt.xlabel("Input V [V]")
plt.ylabel("Output V [V]")
plt.title("Find the linear relationship: a sketch of Output Voltage vs Input Voltage under 100 Hz")
plt.legend()
plt.show()

print(f"CH1 voltage list: {V1}")
print(f"Uncertainty of CH1 voltage list: {uV1}")
print(f"CH2 voltage list: {V2}")
print(f"Uncertainty of CH2 voltage list: {uV2}")

CH1 voltage list: [0.204, 0.306, 0.4, 0.5, 0.598, 0.7, 0.796, 0.892, 0.992, 1.09, 1.19, 1.33, 1.43, 1.52]
Uncertainty of CH1 voltage list: [0.003, 0.003, 0.003, 0.003, 0.003, 0.003, 0.003, 0.003, 0.003, 0.003, 0.003, 0.003, 0.003, 0.003]
CH2 voltage list: [0.96, 1.46, 1.97, 2.5, 2.96, 3.47, 3.97, 4.43, 4.92, 5.44, 5.92, 6.56, 6.92, 7.24]
Uncertainty of CH2 voltage list: [0.020000000000000018, 0.020000000000000018, 0.015000000000000013, 0.014999999999999902, 0.010000000000000009, 0.010000000000000009, 0.010000000000000009, 0.010000000000000231, 0.020000000000000018, 0.020000000000000018, 0.020000000000000018, 0.04, 0.04, 0.04]

de = data_entry2.sheet("finalP-part1amp.csv")

Sheet name: finalP-part1amp.csv

# The script below fits to a sine wave, but can be modified for other functions
# First, we load some python packages
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

################################################
# LIST OF ALL INPUTS
################################################

# fname is assumed to be in a four-column .csv file (comma separated values). The first two rows are
# assumed to be headers, like those produced in our code for packing oscilloscope data.
# The four columns are x-values, x-uncertainties, y-values, y-uncertainties.
# The .csv file must be in the same
# folder as this fit program, otherwise the full file extension must be added
# to fname: e.g. fname = 'folder/subfolder/subsubfolder/file.csv'

fname = "finalP-part1amp.csv"

x_name = "V_input"
x_units = "V"
y_name = "V_output"
y_units = "V"

# The model you will fit to is defined below, in this case a sine wave.
# The parameters in the model are amplitude, freqency, and phase.
# To get a least squares fitting process started, it is necessary to provide good
# initial guesses for the parameters. From your plots of the data so far, you can make good guesses at these parameters.

param_names = ["slope"]


# definition of the fit function
# def fit_function(x, amplitude,tau):

# fit function is a linear model with slope and intercept

def fit_function(x, slope):
    return x*slope

# load the file "fname", defined above# load the file "fname", defined above; 
# if you have file with more than 4 columns, you may need to select different indicies for "usecols" (starting counting at 0)

data = np.loadtxt(fname, delimiter=",", comments="#", usecols=(0, 1, 2, 3), skiprows=2)

################################################
# READ IN DATA COLUMNS
################################################

# Here is where you access the data columns.
# You may need to alter these to choose what is on the y-axis and what is on the x-axis.

x = data[:12, 0]
y = data[:12, 2]
y_sigma = data[:12, 3]

# calculate the best fit slope analytically (1-parameter solution)
m = sum(x*y/y_sigma**2)/sum((x/y_sigma)**2)

# calculate uncertainty of the best fit slope
m_sigma = np.sqrt(1/sum((x/y_sigma)**2))

print ("Slope is", m, "+/-", m_sigma)

###############################################################################
# calculates and prints the chi-squared, degrees of freedon, and weighted chi-squared
###############################################################################

# function that  calculates the chi square value of a fit
def chi_square (param1, x, y, sigma):
#
    return np.sum((y-fit_function(x, param1))**2/sigma**2)
    
# calculate and print chi square as well as the per degree-of-freedom value
chi2 = chi_square(m,x,y,y_sigma)
# degrees of freedom is the number of data points minus the number of parameters
dof = len(x) - 1
print ("\nGoodness of fit - Chi-squared measure:")
print ("degrees of freedom = {}, Chi2/dof = {}\n".format(dof, chi2/dof))


# residual is the difference between the data and model
x_fitfunc = np.linspace(min(x), max(x), 500)
y_fitfunc = fit_function(x_fitfunc, m)
y_fit = fit_function(x, m)
residual = y-y_fit
# creates a histogram of the residuals
hist,bins = np.histogram(residual,bins=30)

fig = plt.figure(figsize=(7,15))

ax1 = fig.add_subplot(311)
ax1.errorbar(x,y,yerr=y_sigma,marker='.',linestyle='',label="measured data")
ax1.plot(x_fitfunc,y_fitfunc,marker="",linestyle="-",linewidth=2,color="r",
         label=" fit")
# add axis labels and title
ax1.set_xlabel('{} [{}]'.format(x_name,x_units))
ax1.set_ylabel('{} [{}]'.format(y_name,y_units))
ax1.set_title('Best fit of 1-Parameter Linear Model')
# set the x and y boundaries of your plot
#plt.xlim(lower_x,upper_x)
#plt.ylim(lower_y,upper_y)
# show a legend. loc='best' places legend where least amount of data is 
# obstructed. 
ax1.legend(loc='best',numpoints=1)


# this code produces a figure with a plot of the residuals as well
# as a histogram of the residuals. 
# fig = plt.figure(figsize=(7,10))
ax2 = fig.add_subplot(312)
ax2.errorbar(x,residual,yerr=y_sigma,marker='.',linestyle='',
             label="residual (y-y_fit)")
ax2.hlines(0,np.min(x),np.max(x),lw=2,alpha=0.8)
ax2.set_xlabel('{} [{}]'.format(x_name,x_units))
ax2.set_ylabel('y-y_fit [{}]'.format(y_units))
ax2.set_title('Residuals for the Best Fit')
ax2.legend(loc='best',numpoints=1)
ax3 = fig.add_subplot(313)
ax3.bar(bins[:-1],hist,width=bins[1]-bins[0])

ax3.set_ylim(0,1.2*np.max(hist))
ax3.set_xlabel('y-y_fit [{}]'.format(y_units))
ax3.set_ylabel('Number of occurences')
ax3.set_title('Histogram of the Residuals')

"""
Modify the following lines to change the name of the file used to store a JPEG of your best fit graphs
"""

# Before showing the plot, you can also save a copy of the figure as a JPEG.
# The order is important here because plt.show clears the plot information after displaying it. 
plt.savefig('FittingResults.jpeg')
plt.show()

Slope is 4.966730980261918 +/- 0.005344811416106966

Goodness of fit - Chi-squared measure:
degrees of freedom = 11, Chi2/dof = 2.369776203085998

# The script below fits to a sine wave, but can be modified for other functions
# First, we load some python packages
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

################################################
# LIST OF ALL INPUTS
################################################

# fname is assumed to be in a four-column .csv file (comma separated values). The first two rows are
# assumed to be headers, like those produced in our code for packing oscilloscope data.
# The four columns are x-values, x-uncertainties, y-values, y-uncertainties.
# The .csv file must be in the same
# folder as this fit program, otherwise the full file extension must be added
# to fname: e.g. fname = 'folder/subfolder/subsubfolder/file.csv'

fname = "finalP-part1amp.csv"

x_name = "V_input"
x_units = "V"
y_name = "V_output"
y_units = "V"

# The model you will fit to is defined below, in this case a sine wave.
# The parameters in the model are amplitude, freqency, and phase.
# To get a least squares fitting process started, it is necessary to provide good
# initial guesses for the parameters. From your plots of the data so far, you can make good guesses at these parameters.

param_names = ["slope", "intercept"]


# definition of the fit function
# def fit_function(x, amplitude,tau):

# fit function is a linear model with slope and intercept

def fit_function(x, slope, intercept):
    return x*slope + intercept

# load the file "fname", defined above; 
# if you have file with more than 4 columns, you may need to select different indicies for "usecols" (starting counting at 0)

data = np.loadtxt(fname, delimiter=",", comments="#", usecols=(0, 1, 2, 3), skiprows=2)

################################################
# READ IN DATA COLUMNS
################################################

# Here is where you access the data columns.
# You may need to alter these to choose what is on the y-axis and what is on the x-axis.

x = data[:12, 0]
y = data[:12, 2]
y_sigma = data[:12, 3]

# calculate the best fit slope - this is the analytical solution for a 2-parameter fit to a striaght line
Delta = sum(1/y_sigma**2)*sum((x/y_sigma)**2) - (sum(x/y_sigma**2)**2)
m = (sum(1/y_sigma**2)*sum(x*y/y_sigma**2) - sum(x/y_sigma**2)*sum(y/y_sigma**2))/Delta

# calculate uncertainty of the best fit slope
m_sigma = np.sqrt(sum(1/y_sigma**2)/Delta)

print ("Slope is", m, "+/-", m_sigma)

# calculate the best fit intercept
b = (sum((x/y_sigma)**2)*sum(y/y_sigma**2) - sum(x/y_sigma**2)*sum(x*y/y_sigma**2))/Delta

# calculate uncertainty of the best fit intercept
b_sigma = np.sqrt(sum((x/y_sigma)**2)/Delta)
                        
print ("Intercept is", b, "+/-", b_sigma)
###############################################################################
# calculates and prints the chi-squared, degrees of freedon, and weighted chi-squared
###############################################################################

# function that  calculates the chi square value of a fit
def chi_square (param1, param2, x, y, sigma):
#
    return np.sum((y-fit_function(x, param1, param2))**2/sigma**2)
    
# calculate and print chi square as well as the per degree-of-freedom value
chi2 = chi_square(m,b,x,y,y_sigma)
# degrees of freedom is the number of data points minus the number of parameters
dof = len(x) - 2
print ("\nGoodness of fit - Chi-squared measure:")
print ("degrees of freedom = {}, Chi2/dof = {}\n".format(dof, chi2/dof))


# residual is the difference between the data and model
x_fitfunc = np.linspace(min(x), max(x), 500)
y_fitfunc = fit_function(x_fitfunc, m ,b)
y_fit = fit_function(x, m,b)
residual = y-y_fit
# creates a histogram of the residuals
hist,bins = np.histogram(residual,bins=30)

fig = plt.figure(figsize=(7,15))

ax1 = fig.add_subplot(311)
ax1.errorbar(x,y,yerr=y_sigma,marker='.',linestyle='',label="measured data")
ax1.plot(x_fitfunc,y_fitfunc,marker="",linestyle="-",linewidth=2,color="r",
         label=" fit")
# add axis labels and title
ax1.set_xlabel('{} [{}]'.format(x_name,x_units))
ax1.set_ylabel('{} [{}]'.format(y_name,y_units))
ax1.set_title('Best fit of 2-Parameter Linear Model')
# set the x and y boundaries of your plot
#plt.xlim(lower_x,upper_x)
#plt.ylim(lower_y,upper_y)
# show a legend. loc='best' places legend where least amount of data is 
# obstructed. 
ax1.legend(loc='best',numpoints=1)


# this code produces a figure with a plot of the residuals as well
# as a histogram of the residuals. 
# fig = plt.figure(figsize=(7,10))
ax2 = fig.add_subplot(312)
ax2.errorbar(x,residual,yerr=y_sigma,marker='.',linestyle='',
             label="residual (y-y_fit)")
ax2.hlines(0,np.min(x),np.max(x),lw=2,alpha=0.8)
ax2.set_xlabel('{} [{}]'.format(x_name,x_units))
ax2.set_ylabel('y-y_fit [{}]'.format(y_units))
ax2.set_title('Residuals for the Best Fit')
ax2.legend(loc='best',numpoints=1)
ax3 = fig.add_subplot(313)
ax3.bar(bins[:-1],hist,width=bins[1]-bins[0])

ax3.set_ylim(0,1.2*np.max(hist))
ax3.set_xlabel('y-y_fit [{}]'.format(y_units))
ax3.set_ylabel('Number of occurences')
ax3.set_title('Histogram of the Residuals')

"""
Modify the following lines to change the name of the file used to store a JPEG of your best fit graphs
"""

# Before showing the plot, you can also save a copy of the figure as a JPEG.
# The order is important here because plt.show clears the plot information after displaying it. 
plt.savefig('FittingResults.jpeg')
plt.show()

Slope is 5.0156414179820885 +/- 0.016716139153860488
Intercept is -0.038783139148621706 +/- 0.012559117049106876

Goodness of fit - Chi-squared measure:
degrees of freedom = 10, Chi2/dof = 1.653150621078661

# 1kHz
uVin1 = (196e-3 - 193e-3) / 2  #V
uVout1 = 50e-3 / 2  #V
G1 = 960e-3 / 196e-3  # Vout/Vin
uG1 = G1 * np.sqrt((uVin1 / 196e-3)**2 + (uVout1 / 960e-3)**2)

# 3kHz
uVin2 = (197e-3 - 195e-3) / 2  #V
uVout2 = (980e-3 - 960e-3) / 2  #V
G2 = 980e-3 / 197e-3  # Vout/Vin
uG2 = G2 * np.sqrt((uVin2 / 197e-3)**2 + (uVout2 / 980e-3)**2)

# 6kHz
uVin3 = (197e-3 - 195e-3) / 2  #V
uVout3 = (960e-3 - 940e-3) / 2  #V
G3 = 960e-3 / 197e-3  # Vout/Vin
uG3 = G3 * np.sqrt((uVin3 / 197e-3)**2 + (uVout3 / 960e-3)**2)

# 10kHz
uVin4 = (195e-3 - 193e-3) / 2  #V
uVout4 = (960e-3 - 940e-3) / 2  #V
G4 = 960e-3 / 195e-3  # Vout/Vin
uG4 = G4 * np.sqrt((uVin4 / 195e-3)**2 + (uVout4 / 960e-3)**2)

# 14kHz
uVin5 = (195e-3 - 193e-3) / 2  #V
uVout5 = (960e-3 - 940e-3) / 2  #V
G5 = G4  # same data of Vin and Vout
uG5 = uG4

# 20kHz
uVin6 = (193e-3 - 191e-3) / 2  #V
uVout6 = (940e-3 - 920e-3) / 2  #V
G6 = 940e-3 / 193e-3  # Vout/Vin
uG6 = G6 * np.sqrt((uVin6 / 193e-3)**2 + (uVout6 / 940e-3)**2)

# 30kHz
uVin7 = (192e-3 - 191e-3) / 2  #V
uVout7 = (920e-3 - 900e-3) / 2  #V
G7 = 920e-3 / 191e-3  # Vout/Vin
uG7 = G7 * np.sqrt((uVin7 / 191e-3)**2 + (uVout7 / 920e-3)**2)

# 35kHz
uVin8 = (191e-3 - 189e-3) / 2  #V
uVout8 = (940e-3 - 880e-3) / 2  #V
G8 = 920e-3 / 188e-3  # Vout/Vin
uG8 = G8 * np.sqrt((uVin8 / 188e-3)**2 + (uVout8 / 920e-3)**2)

# 40kHz
uVin9 = (188e-3 - 186e-3) / 2  #V
uVout9 = (920e-3 - 880e-3) / 2  #V
G9 = 900e-3 / 186e-3  # Vout/Vin
uG9 = G9 * np.sqrt((uVin9 / 186e-3)**2 + (uVout9 / 900e-3)**2)

# 45kHz
uVin10 = (186e-3 - 184e-3) / 2  # V
uVout10 = (930e-3 - 870e-3) / 2  # V
G10 = 900e-3 / 184e-3  # Vout/Vin
uG10 = G10 * np.sqrt((uVin10 / 184e-3)**2 + (uVout10 / 900e-3)**2)

# 100kHz
uVin11 = (155e-3 - 153e-3) / 2  #V
uVout11 = (770e-3 - 710e-3) / 2  #V
G11 = 770e-3 / 154e-3  # Vout/Vin
uG11 = G11 * np.sqrt((uVin11 / 154e-3)**2 + (uVout11 / 770e-3)**2)

# 150 -> same for 200, 250 kHz
uVin12 = 1e-3 / 2  #V
uVout12 = 0.1 / 2  #V
G12 = 610e-3 / 126e-3  # Vout/Vin
uG12 = G12 * np.sqrt((uVin12 / 126e-3)**2 + (uVout12 / 610e-3)**2)

# 200kHz
uVin13 = uVin12  #V
uVout13 = uVout12  #V
G13 = 510e-3 / 106e-3  # Vout/Vin
uG13 = G13 * np.sqrt((uVin13 / 106e-3)**2 + (uVout13 / 510e-3)**2)

# 250kHz
uVin14 = uVin13  #V
uVout14 = uVout13  #V
G14 = 410e-3 / 86e-3  # Vout/Vin
uG14 = G14 * np.sqrt((uVin14 / 86e-3)**2 + (uVout14 / 410e-3)**2)

# 24kHz
uVin15 = uVin6 #V
uVout15 = uVout6
G15 = 940e-3/193e-3 # Vout/Vin
uG15 = G15 * np.sqrt((uVin15 / 193e-3)**2 + (uVout15 / 940e-3)**2)

uV1 = [uVin1, uVin2, uVin3, uVin4, uVin5, uVin6, uVin15, uVin7, uVin8, uVin9, uVin10, uVin11, uVin12, uVin13, uVin14]
uV2 = [uVout1, uVout2, uVout3, uVout4, uVout5, uVout6, uVout15, uVout7, uVout8, uVout9, uVout10, uVout11, uVout12, uVout13, uVout14]
print(f"Uncertainty of CH1 voltage list: {uV1}")
print(f"Uncertainty of CH2 voltage list: {uV2}")

G = [G1, G2, G3, G4, G5, G6, G15, G7, G8, G9, G10, G11, G12, G13, G14]
uG = [uG1, uG2, uG3, uG4, uG5, uG6, uG15, uG7, uG8,uG9, uG10, uG11, uG12, uG13, uG14]

Uncertainty of CH1 voltage list: [0.0015000000000000013, 0.0010000000000000009, 0.0010000000000000009, 0.0010000000000000009, 0.0010000000000000009, 0.0010000000000000009, 0.0010000000000000009, 0.0005000000000000004, 0.0010000000000000009, 0.0010000000000000009, 0.0010000000000000009, 0.0010000000000000009, 0.0005, 0.0005, 0.0005]
Uncertainty of CH2 voltage list: [0.025, 0.010000000000000009, 0.010000000000000009, 0.010000000000000009, 0.010000000000000009, 0.009999999999999953, 0.009999999999999953, 0.010000000000000009, 0.02999999999999997, 0.020000000000000018, 0.030000000000000027, 0.030000000000000027, 0.05, 0.05, 0.05]

de = data_entry2.sheet("finalP-part1freq.csv")

Sheet name: finalP-part1freq.csv

de = data_entry2.sheet("finalP-part1freqnew.csv")

Sheet name: finalP-part1freqnew.csv

plt.errorbar(FreqVec,VoutVec, yerr=dVoutVec, marker=".", linestyle="",label="Measured data")
plt.xlabel("Frequency(input) [kHz]")
plt.ylabel("Vout [V]")
plt.title("Find the relationship: a sketch of the Amplitude of output siganl versus the Frequency of Input signal")
plt.legend()
plt.show()

C3 = 0.425 # nF #9.999 nF, so resolution is 0.001 nF, accuracy is 5%+20*unit
dC3 = 0.425*0.05+20*(0.001)

R4 = 100.3e3 # Ohms #999.9 kOhms, so resolution is 0.1kOhms, accuracy is 0.5%+3*unit
dR4 = 100.3e3*0.005+3*(0.1*1000)

# You may notice that i use the value of some uncertainty directly, because i want to hold all values in one counting method
print(f"DMM: 100 fixed kOhms(R4) resistor is {R4:.0f} +/- 800 Ohms")
print(f"DMM: 0.5nF capacitor(C3) is {C3:.2f} +/- {dC3:.1g} nF")

DMM: 100 fixed kOhms(R4) resistor is 100300 +/- 800 Ohms
DMM: 0.5nF capacitor(C3) is 0.42 +/- 0.04 nF

# INPUT SIGNAL measured by cursors
amp_in = 90e-3 #V
damp_in = (110e-3-60e-3)/2
freq_in = 692 #kHz
dfreq_in = (716-666)/2

print(f"the amplitude of the input signal is {amp_in} +/- {damp_in:.1g} V")
print(f"the frequency of the input signal is 690 +/- 20 kHz")

the amplitude of the input signal is 0.09 +/- 0.03 V
the frequency of the input signal is 690 +/- 20 kHz

# OUTPUT SIGNAL: compare the amp and freq on the right side of the diode
amp_out = 600e-3 #V
damp_out = (812e-3-560e-3)/2
freq_out = 6.848 #kHz
dfreq_out = (6.972-6.734)/2
Voffset_out = 460e-3 #V
dVoffset_out = (500e-3-400e-3)/2

print(f"the amplitude of the output signal is {amp_out} +/- {damp_out:.1g} V")
print(f"the frequency of the output signal is {freq_out:.1f} +/- {dfreq_out:.1g} kHz")
print(f"the offset of the output signal is {Voffset_out} +/- {dVoffset_out:.1g} V")

the amplitude of the output signal is 0.6 +/- 0.1 V
the frequency of the output signal is 6.8 +/- 0.1 kHz
the offset of the output signal is 0.46 +/- 0.05 V

# OUTPUT SIGNAL
amp_out = 3.06 #V
damp_out = (3.16-2.80)/2
freq_out = 6.900 #kHz
dfreq_out = (6.922-6.867)/2
Voffset_out = 0 #V

print(f"the amplitude of the output signal is {amp_out:.1f} +/- {damp_out:.1g} V")
print(f"the frequency of the output signal is {freq_out:.2f} +/- {dfreq_out:.1g} kHz")
print(f"the offset of the output signal is 0 V")

the amplitude of the output signal is 3.1 +/- 0.2 V
the frequency of the output signal is 6.90 +/- 0.03 kHz
the offset of the output signal is 0 V

L1 = 450 #uH
dL1 = (450 - 449)/2 #uH

R5 = 98.79 # Ohms #99.99 Ohms, so resolution is 0.01 Ohms, accuracy is 1%+3*unit
dR5 = 98.79*0.01+3*(0.01)

print(f"The new 100 Ohms resistance is {R5:.0f} +/- {dR5:.1g} Ohms ")
print(f"The inductance is {L1:.1f} +/- {dL1:.1g} uH")

The new 100 Ohms resistance is 99 +/- 1 Ohms 
The inductance is 450.0 +/- 0.5 uH

# Freq_input=690 kHz
f = 690 # kHz

V_min = 6.2 # V
V_max = 6.4 # V
V = 6.2 # V

# Calculate the uncertainty, using Gaussian distribution and 68% certain interval
uf =  0.1/2
uV = (V_max - V_min)/2 
print(f"f1 = {f:.2f} +/- {uf:.1g}, and V1 = {V:.2f} +/- {uV:.1g}.")

# Freq_input=590 kHz
f2 = 590 # kHz

V_min2 = 1.6 # V
V_max2 = 1.84 # V
V2 = 1.78 # V

# Calculate the uncertainty
uf2 = 0.1/2
uV2 = (V_max2 - V_min2)/2 
print(f"f2 = {f2:.2f} +/- {uf2:.1g}, and V2 = {V2:.2f} +/- {uV2:.1g}.")

# Freq_input=600 kHz
f3 = 600 # kHz

V_min3 = 1.76 # V
V_max3 = 1.96 # V
V3 = 1.94 # V

# Calculate the uncertainty
uf3 = 0.1/2
uV3 = (V_max3 - V_min3)/2 
print(f"f3 = {f3:.2f} +/- {uf3:.1g}, and V3 = {V3:.2f} +/- {uV3:.1g}.")

# Freq_input=610 kHz
f4 = 610 # kHz

V_min4 = 1.86 # V
V_max4 = 2.22 # V
V4 = 2.12 # V

# Calculate the uncertainty
uf4 = 0.1/2
uV4 = (V_max4 - V_min4)/2 
print(f"f4 = {f4:.2f} +/- {uf4:.1g}, and V4 = {V4:.2f} +/- {uV4:.1g}.")

# Freq_input=620 kHz
f5 = 620 # kHz

V_min5 = 2.18 # V
V_max5 = 2.56 # V
V5 = 2.46  # V

# Calculate the uncertainty
uf5 = 0.1/2
uV5 = (V_max5 - V_min5)/2 
print(f"f5 = {f5:.2f} +/- {uf5:.1g}, and V5 = {V5:.2f} +/- {uV5:.1g}.")

# Freq_input=630 kHz
f6 = 630 # kHz

V_min6 = 2.76 # V
V_max6 =2.92 # V
V6 = 2.92 # V

# Calculate the uncertainty
uf6 = 0.1/2
uV6 = (V_max6 - V_min6)/2 
print(f"f6 = {f6:.2f} +/- {uf6:.1g}, and V6 = {V6:.2f} +/- {uV6:.1g}.")

# Freq_input=640 kHz
f7 = 641.0 # kHz

V_min7 = 3.00 # V
V_max7 = 3.16  # V
V7 = 3.16 # V

# Calculate the uncertainty
uf7 = 0.1/2
uV7 = (V_max7 - V_min7)/2 
print(f"f7 = {f7:.2f} +/- {uf7:.1g}, and V7 = {V7:.2f} +/- {uV7:.1g}.")

# Freq_input=650 kHz
f8 = 655 # kHz

V_min8 = 3.56 # V
V_max8 = 3.72 # V
V8 = 3.72  # V

# Calculate the uncertainty
uf8 = 0.1/2
uV8 = (V_max8 - V_min8)/2 
print(f"f8 = {f8:.2f} +/- {uf8:.1g}, and V8 = {V8:.2f} +/- {uV8:.1g}.")

# Freq_input=700 kHz
f9 = 700.0 # kHz

V_min9 = 6.36 # V
V_max9 = 6.44 # V
V9 = 6.44  # V

# Calculate the uncertainty
uf9 = 0.1/2
uV9 = (V_max9 - V_min9)/2 
print(f"f9 = {f9:.2f} +/- {uf9:.1g}, and V9 = {V9:.2f} +/- {uV9:.1g}.")

# Freq_input=710 kHz
f10 = 710 # kHz

V_min10 = 6.36 # V
V_max10 = 6.60 # V
V10 = 6.46 # V

# Calculate the uncertainty
uf10 = 0.1/2
uV10 = (V_max10 - V_min10)/2 
print(f"f10 = {f10:.2f} +/- {uf10:.1g}, and V10 = {V10:.2f} +/- {uV10:.1g}.")

# Freq_input= 680kHz
f11 = 680 # kHz

V_min11 = 5.48 # V
V_max11 = 5.72 # V
V11 = 5.56 # V

# Calculate the uncertainty
uf11 = 0.1/2
uV11 = (V_max11 - V_min11)/2 
print(f"f11 = {f11:.1f} +/- {uf11:.1g}, and V11 = {V11:.2f} +/- {uV11:.1g}.")

# Freq_input=685 kHz
f12 = 685 # kHz

V_min12 = 5.72 # V
V_max12 = 6.00 # V
V12 = 5.80 # V

# Calculate the uncertainty
uf12 = 0.1/2
uV12 = (V_max12 - V_min12)/2 
print(f"f12 = {f12:.2f} +/- {uf12:.1g}, and V12 = {V12:.2f} +/- {uV12:.1g}.")

# Freq_input=670 kHz
f13 = 670 # kHz

V_min13 = 4.84 # V
V_max13 = 5.00 # V
V13 = 4.92 # V

# Calculate the uncertainty
uf13 = 0.1/2
uV13 = (V_max13 - V_min13)/2 
print(f"f13 = {f13:.1f} +/- {uf13:.1g}, and V13 = {V13:.2f} +/- {uV13:.1g}.")

# Freq_input=720 kHz
f14 = 720 # kHz

V_min14 = 5.92 # V
V_max14 = 6.04 # V
V14 = 6.04 # V

# Calculate the uncertainty
uf14 = 0.1/2
uV14 = (V_max14 - V_min14)/2 
print(f"f14 = {f14:.2f} +/- {uf14:.1g}, and V14 = {V14:.2f} +/- {uV14:.1g}.")

# Freq_input=750 kHz
f15 = 750 # kHz

V_min15 = 4.12 # V
V_max15 = 4.32 # V
V15 = 4.20 # V

# Calculate the uncertainty
uf15 = 0.1/2
uV15 = (V_max15 - V_min15)/2 
print(f"f15 = {f15:.1f} +/- {uf15:.1g}, and V15 = {V15:.2f} +/- {uV15:.1g}.")

# Freq_input=780 kHz
f16 = 780 # kHz

V_min16 = 2.56 # V
V_max16 = 2.92 # V
V16 = 2.84 # V

# Calculate the uncertainty
uf16 = 0.1/2
uV16 = (V_max16 - V_min16)/2 
print(f"f16 = {f16:.2f} +/- {uf16:.1g}, and V16 = {V16:.2f} +/- {uV16:.1g}.")

# Freq_input=730 kHz
f17 = 730 # kHz

V_min17 = 5.40 # V
V_max17 = 5.56 # V
V17 = 5.48 # V

# Calculate the uncertainty
uf17 = 0.1/2
uV17 = (V_max17 - V_min17)/2 
print(f"f17 = {f17:.2f} +/- {uf17:.1g}, and V17 = {V17:.2f} +/- {uV17:.1g}.")

# Freq_input= 740kHz
f18 = 740 # kHz

V_min18 = 4.76 # V
V_max18 = 5.04 # V
V18 = 4.84 # V

# Calculate the uncertainty
uf18 = 0.1/2
uV18 = (V_max18 - V_min18)/2 
print(f"f18 = {f18:.2f} +/- {uf18:.1g}, and V18 = {V18:.2f} +/- {uV18:.1g}.")

# Freq_input=760 kHz
f19 = 760 # kHz

V_min19 = 3.64 # V
V_max19 = 3.80 # V
V19 = 3.72 # V

# Calculate the uncertainty
uf19 = 0.1/2
uV19 = (V_max19 - V_min19)/2 
print(f"f19 = {f19:.1f} +/- {uf19:.1g}, and V19 = {V19:.2f} +/- {uV19:.1g}.")

# Freq_input=800 kHz
f20 = 800 # kHz

V_min20 = 2.26 # V
V_max20 = 2.36 # V
V20 = 2.28 # V

# Calculate the uncertainty
uf20 = 0.1/2
uV20 = (V_max20 - V_min20)/2 
print(f"f20 = {f20:.1f} +/- {uf20:.1g}, and V20 = {V20:.2f} +/- {uV20:.1g}.")

# Freq_input=850 kHz
f21 = 850 # kHz

V_min21 = 1.48 # V
V_max21 = 1.60 # V
V21 = 1.58 # V

# Calculate the uncertainty
uf21 = 0.1/2
uV21 = (V_max21 - V_min21)/2 
print(f"f21 = {f21:.1f} +/- {uf21:.1g}, and V21 = {V21:.2f} +/- {uV21:.1g}.")

# Freq_input=550 kHz
f22 = 550 # kHz

V_min22 = 1.08 # V
V_max22 = 1.32 # V
V22 = 1.26 # V

# Calculate the uncertainty
uf22 = 0.1/2
uV22 = (V_max22 - V_min22)/2 
print(f"f22 = {f22:.1f} +/- {uf22:.1g}, and V22 = {V22:.2f} +/- {uV22:.1g}.")

V = [V, V2, V3, V4, V5, V6, V7, V8, V9, V10, V11, V12, V13, V14, V15, V16, V17, V18, V19, V20, V21, V22]
f = [f, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12, f13, f14, f15, f16, f17, f18, f19, f20, f21, f22]
uV = [uV, uV2, uV3, uV4, uV5, uV6, uV7, uV8, uV9, uV10, uV11, uV12, uV13, uV14, uV15, uV16, uV17, uV18, uV19, uV20, uV21, uV22]
plt.errorbar(f,V,yerr=uV,marker='.',linestyle='',label="LRC data")
plt.title("Find the relationship: a sketch of the amplitude of output signal versus the frequency of input signal")
plt.xlabel("Frequency [kHz]")
plt.ylabel("Amplitude [V]")
plt.legend()
plt.show()

f1 = 690.00 +/- 0.05, and V1 = 6.20 +/- 0.1.
f2 = 590.00 +/- 0.05, and V2 = 1.78 +/- 0.1.
f3 = 600.00 +/- 0.05, and V3 = 1.94 +/- 0.1.
f4 = 610.00 +/- 0.05, and V4 = 2.12 +/- 0.2.
f5 = 620.00 +/- 0.05, and V5 = 2.46 +/- 0.2.
f6 = 630.00 +/- 0.05, and V6 = 2.92 +/- 0.08.
f7 = 641.00 +/- 0.05, and V7 = 3.16 +/- 0.08.
f8 = 655.00 +/- 0.05, and V8 = 3.72 +/- 0.08.
f9 = 700.00 +/- 0.05, and V9 = 6.44 +/- 0.04.
f10 = 710.00 +/- 0.05, and V10 = 6.46 +/- 0.1.
f11 = 680.0 +/- 0.05, and V11 = 5.56 +/- 0.1.
f12 = 685.00 +/- 0.05, and V12 = 5.80 +/- 0.1.
f13 = 670.0 +/- 0.05, and V13 = 4.92 +/- 0.08.
f14 = 720.00 +/- 0.05, and V14 = 6.04 +/- 0.06.
f15 = 750.0 +/- 0.05, and V15 = 4.20 +/- 0.1.
f16 = 780.00 +/- 0.05, and V16 = 2.84 +/- 0.2.
f17 = 730.00 +/- 0.05, and V17 = 5.48 +/- 0.08.
f18 = 740.00 +/- 0.05, and V18 = 4.84 +/- 0.1.
f19 = 760.0 +/- 0.05, and V19 = 3.72 +/- 0.08.
f20 = 800.0 +/- 0.05, and V20 = 2.28 +/- 0.05.
f21 = 850.0 +/- 0.05, and V21 = 1.58 +/- 0.06.
f22 = 550.0 +/- 0.05, and V22 = 1.26 +/- 0.1.

de = data_entry2.sheet("finalP-part3freq.csv")

Sheet name: finalP-part3freq.csv

# Load python packages
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit

###############################################################################
# DEFINED FITTING FUNCTIONS
###############################################################################

def sine_func(x, amplitude, freq, phase):
    return amplitude * np.sin(2.0 * np.pi * freq * x + phase)

def offset_sine_func(x, amplitude, freq, phase, offset):
    return (amplitude * np.sin(2.0 * np.pi * freq * x + phase)) + offset

def exponential_func(x, amplitude, tau, voffset):
    return amplitude * np.exp(x/(-1.0*tau)) + voffset

def ringdown_function(x, amplitude, tau, resonantf, phase):
    return amplitude * np.exp(-x/tau) * np.cos(2.0*np.pi * resonantf * x + phase)

def linear_func(x, slope, intercept):
    return slope * x + intercept

def RCresp_func(x, tau, V0, voffset):
    return V0/np.sqrt(1+(2*np.pi*tau*x)**2) + voffset

def LRCresp_func(x, f0, Vin, gamma, foffset):
    return (Vin/(np.sqrt(1 + (2 * np.pi / (gamma * (x + foffset)))**2 * ((x + foffset)**2 - f0**2)**2)))

def LRCresp_func2(x, f0, Vin, gamma, Voffset):
    return (Vin/(np.sqrt(1 + (2 * np.pi / (gamma * (x)))**2 * ((x)**2 - f0**2)**2))) + Voffset

def LRCresp_func3(x, f0, Vin, gamma):
    return (Vin/(np.sqrt(1 + (2 * np.pi / (gamma * (x)))**2 * ((x)**2 - f0**2)**2)))

###############################################################################
# LIST OF ALL INPUTS
###############################################################################

# Name of the data file
fname = "finalP-part3freq.csv"

# Names and units of data columns from fname
x_name = "frequency"
x_units = "kHz"
y_name = "Voltage"
y_units = "V"

# Modify to change the fitting function, parameter names and to set initial parameter guesses
fit_function = LRCresp_func2
param_names = ("f0", "Vin", "gamma", "Voffset")
guesses = (744,4.08,50, 0)

# Flags for optional features
show_covariance_matrix = False
set_xy_boundaries = False
lower_x = -0.01 # these values ignored if set_xy_boundaries = False
upper_x = 0.01
lower_y = -1
upper_y = 1

###############################################################################
# LOAD DATA
###############################################################################

# load the file fname and skip the first 'skiprows' rows
data = np.loadtxt(fname, delimiter=",", comments="#", usecols=(0, 1, 2, 3), skiprows=2)

# Assign the data file columns to variables for later use
x = data[:, 0]
y = data[:, 2]
y_sigma = data[:, 3]

###############################################################################
# INITIAL PLOT OF THE DATA
###############################################################################

# Define 500 points spanning the range of the x-data; for plotting smooth curves
xtheory = np.linspace(min(x), max(x), 500)

# Compare the guessed curve to the data for visual reference
y_guess = fit_function(xtheory, *guesses)
plt.errorbar(x, y, yerr=y_sigma, marker=".", linestyle="", label="Measured data")
plt.plot(
    xtheory,
    y_guess,
    marker="",
    linestyle="-",
    linewidth=1,
    color="g",
    label="Initial parameter guesses",
)
plt.xlabel(f"{x_name} [{x_units}]")
plt.ylabel(f"{y_name} [{y_units}]")
plt.title(r"Comparison between the data and the intial parameter guesses")
plt.legend(loc="best", numpoints=1)
plt.show()

# calculate the value of the model at each of the x-values of the data set
y_fit = fit_function(x, *guesses)

# Residuals are the difference between the data and theory
residual = y - y_fit

# Plot the residuals
plt.errorbar(x, residual, yerr=y_sigma, marker=".", linestyle="", label="residuals")
plt.xlabel(f"{x_name} [{x_units}]")
plt.ylabel(f"Residual y-y_fit [{y_units}]")
plt.title("Residuals using initial parameter guesses")
plt.show()

###############################################################################
# PERFORM THE FIT AND PRINT RESULTS
###############################################################################

# Use curve_fit to perform the fit
# fit_function: defined above to choose a specific fitting function 
# fit_params: holds the resulting fit parameters
# fit_cov: the covariance matrix between all the parameters
#          (used to extract fitting parameter uncertanties)
# maxfev=10**5: maximum number of fitting procedure iterations before giving up
# absolute_sigma=True: uncertainties are treated as absolute (not relative)
fit_params, fit_cov = curve_fit(
    fit_function, x, y, sigma=y_sigma, 
    p0=guesses,absolute_sigma=True, maxfev=10**5)

# Define the function that calculates chi-squared
def chi_square(fit_parameters, x, y, sigma):
    dof = len(x) - len(fit_params)
    return np.sum((y - fit_function(x, *fit_parameters)) ** 2 / sigma**2)/dof

# Calculate and print reduced chi-squared
chi2 = chi_square(fit_params, x, y, y_sigma)
print("Chi-squared = ", chi2)

# Calculate the uncertainties in the fit parameters
fit_params_error = np.sqrt(np.diag(fit_cov))

# Print the fit parameters with uncertianties
print("\nFit parameters:")
for i in range(len(fit_params)):
    print(f"   {param_names[i]} = {fit_params[i]:.3e} ± {fit_params_error[i]:.3e}")
print("\n")

# (Optional) Print the covariance between all variables
if show_covariance_matrix:
    print("Covariance between fit parameters:")
    for i, fit_covariance in enumerate(fit_cov):
        for j in range(i+1,len(fit_covariance)):
            print(f"   {param_names[i]} and {param_names[j]}: {fit_cov[i,j]:.3e}")
    print("\n")

# residual is the difference between the data and model
x_fitfunc = np.linspace(min(x), max(x), 500)
y_fitfunc = fit_function(x_fitfunc, *fit_params)
y_fit = fit_function(x, *fit_params)
residual = y-y_fit

###############################################################################
# PRODUCE A MULTIPANEL PLOT, WITH SCATTER PLOT, RESIDUALS AND RESIDUALS HISTOGRAM
###############################################################################

# The size of the canvas
fig = plt.figure(figsize=(7,15))

# The scatter plot
ax1 = fig.add_subplot(311)
ax1.errorbar(x,y,yerr=y_sigma,marker='.',linestyle='',label="Measured data")
ax1.plot(x_fitfunc, y_fitfunc, marker="", linestyle="-", linewidth=2,color="r", label="Fit")
ax1.set_xlabel(f"{x_name} [{x_units}]")
ax1.set_ylabel(f"{y_name} [{y_units}]")
ax1.set_title('Best Fit of Function to Data')

# (Optional) set the x and y boundaries of your plot
if set_xy_boundaries:
    plt.xlim(lower_x,upper_x)
    plt.ylim(lower_y,upper_y)
# Show the legend. loc='best' places it where the date are least obstructed
ax1.legend(loc='best',numpoints=1)

# The residuals plot
ax2 = fig.add_subplot(312)
ax2.errorbar(x, residual, yerr=y_sigma,marker='.', linestyle='', label="Residual (y-y_fit)")
ax2.hlines(0,np.min(x),np.max(x),lw=2,alpha=0.8)
ax2.set_xlabel(f"{x_name} [{x_units}]")
ax2.set_ylabel(f"y-y_fit [{y_units}]")
ax2.set_title('Residuals for the Best Fit')
ax2.legend(loc='best',numpoints=1)

# Histogram of the residuals
ax3 = fig.add_subplot(313)
hist,bins = np.histogram(residual,bins=30)
ax3.bar(bins[:-1],hist,width=bins[1]-bins[0])
ax3.set_ylim(0,1.2*np.max(hist))
ax3.set_xlabel(f"y-y_fit [{y_units}]")
ax3.set_ylabel('Number of occurences')
ax3.set_title('Histogram of the Residuals')

# Save a copy of the figure as a png 
plt.savefig('FittingResults.png')

# Show the plot
plt.show()

Chi-squared =  1.1891120313124306

Fit parameters:
   f0 = 7.037e+02 ± 4.116e-01
   Vin = 6.913e+00 ± 9.241e-02
   gamma = 5.025e+02 ± 1.332e+01
   Voffset = -4.024e-01 ± 9.828e-02

#2024 Final Project: AM Radio Receiver¶

Background (from instruction and tutorial)¶

Experiment goal¶

The complete circuit¶

During Lab: Part 1: Build and Characterize the Amplifier¶

Experiment equipment¶

Experiment Procedure and Plan¶

The circuit of part 1 we actually built¶

Some Observation (under analysis of part 1)¶

Record data from DMM measured¶

Record data of $V_{output}$ versus $V_{input}$¶

Comment and Analyze fit for $V_{output}$ versus $V_{input}$¶

Record data of adjust frequency of input signal¶

Important Update here: (Date of update: 2024.11.29)¶

Analyze the relationship between the frequency of the input signal and the gain¶

Record data of adjust Voffset¶

During Lab: Part 2: Build and Characterize the Demodulation¶

Experiment equipment¶

Experiment Procedure and Plan¶

The circuit of part 2 we actually built¶

Record data from DMM measured¶

Record data of Amplitude, offset and frequency of the output signal on CH2¶

Analysis of part 2¶

During Lab: Part 3: Build and Characterize the Tank circuit¶

Experiment equipment¶

Experiment Procedure and Plan¶

The circuit of part 3 we actually built¶

Observation¶

Record data from DMM measured¶

Record data of the amplitude of Vout versus the frequency of Vin¶

Comment and Analyze the fit of part 3¶

During Lab: Final Part: Testing the AM receiver¶

Experiment equipment¶

Experiment Procedure and Plan¶

The circuit we finally built¶

Observation and Analysis¶

Summary: Issues and Improvements¶